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LAGRANGIAN SIMULATION OF TAYLOR-COUETTE FLOW 


I. Introduction 

We report here on the development of the hydrodynamics code SPLASH 
which is designed for the Lagrangian simulation of transient rotational flow 
phenomena. The code solves the incompressible, inviscid fluid equations in 
an axisymmetric, cylindrical coordinate system. This is a 2 1/2 dimensional 
model with two spatial coordinates (r,z) and three velocity coordinates 
(u,w,v). All variables are independent of the angle 0. 

The equations of motion are finite-differenced on a general-connectivity 
triangular mesh. A triangular mesh is the natural choice for flows in 
complicated geometries or flows with free surfaces or interfaces. The 
general-connectivity mesh allows local mesh restructuring whenever the grid 
distorts sufficiently to affect numerical accuracy and convergence. The 
SPLASH code is a direct extension of the hydrodynamics code SPLISH^ which 
solves the incompressible, inviscid hydrodynamic equations on a triangular 
mesh in Cartesian geometry. 

This model is applicable to a host of problems such as rotating columns 

of fluid including imploding liner systems with axially displaced annular or 
2 

end plate pistons , axisymmetric jets, laser ablation of spherical shells 

and droplet combustion. Here we apply the model to the study of Couette flow 
3 

and Taylor vortex formation between two rotating coaxial cylinders. This 

problem was chosen as a test case because the linear theory is straight- 

4 5 

forward and well-developed ’ and there is a myriad of experimental 
6 7 8 

results ’ ’ available for comparison. We have been unable to find any 
numerical work in the literature which models the time-evolution of rota¬ 
tional flows to serve as a comparison. 


Manuscript submitted May 14, 1981. 
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In Che next section we discuss several aspects of the triangular 
gridding techniques. The equations of motion are developed in Section Ill 
with the finite-difference algorithms presented in Section IV. In Section V 
we discuss the Couette flow problem and the numerical results are presented 
in Section VI. The summary and conclusions make up Section VII. 


II. Triangular Mesh 

The set of equations governing incompressible, inviscid flow in a 
cylindrical coordinate system will be approximated by finite differences 
on a triangular mesh. The variables in these equations will be represented 
as triangle or vertex quantities on this mesh. This differencing procedure 
is somewhat complex. We will illustrate it by discussing some important 
basic concepts. 

The basic computational cell is the shaded region shown in Fig. 1. It 

is formed by joining the side bisectors of the triangles surrounding the 

general vertex. Each triangle surrounding a central vertex contributes 

1/3 of its area to the area of the basic cell. 

We now illustrate how a gradient is represented in this model. If 

vertex quantities are linear functions of position, then, given the 

function g (defined on vertex m), the function g at any other point, n, 
m 

say, can be written without approximation, as 


8 n = 


m 


+ R 
— n 




(II-l) 


Here R is the vector from the location of g to the chosen point. 

— n m 

Now consider the triangle j defined by two side vectors , _S 
with the vertex-defined quantities g, g^ and The index j indicates 

triangle quantities, the index i vertex quantities (see Fig. 2). Following 
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Fig. 1 — Detail of the triangular grid elements. A triangle is made up of three 
directed line segments, and a vertex represents the shaded region defined by 
the center of mass of each surrounding triangle and the midpoint of each side. 
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Ref. 9, the gradient of g, uniquely defined on the triangle and constant 
throughout it is, 


Vg. - (g r g) s ±+1 


( g i+ l-g) i. 


■ I n 

i=*l 


Ik. 

J 


(II-2) 


where S_ is the side vector rotated clockwise by tt/2 radians, and is short¬ 
hand for a cross product with a basis vector. Here n is a unit vector normal 
to the computational plane and is the area of the triangle. 

2Aj = (jr i x (£.“ JL i+ ^) • n. _r^ is the coordinate location of the 
vertex i. The conclusion is that gradients may be naturally represented on 
triangles, and easily calculated on them in the linear approximation. 

The integral operator consists of a piecewise triangle summation 
about a basic cell giving rise to a vertex integral which is likewise exact 
for the linear approximation. Expressions for the divergence, curl and 
Laplacian are presented in detail in Ref. 1. The finite differencing of 
the nonlinear diffusion equation (Eq. A-5) is discussed in Ref. 10. Although 
the basic difference and integral operators are linear, the resulting 
weighting to a central vertex yields second order accurate approximations 
in the same way that central differences are second order accurate for one- 
dimension. 

Although the control volume approach assures that the equations are 
solved conservatively, large numerical errors may arise in a Lagrangian code 
due to severe grid distortion. If portions of the grid become stretched, 
gradients will be calculated which involve vertices far removed from one 
another. The convergence of the iterations would be slow and truncation 
errors would build rapidly as the triangle sides lengthen. 







This difficulty is avoided by forcing the mesh to restructure. Mesh 
restructuring may involve interchanging the diagonals of a quadrilateral 
formed by two adjacent triangles or adding/deleting a vertex on a triangle 
side or in the interior of a triangle. Restructuring is performed to pre¬ 
serve the diagonal dominance of the Poisson equation as a specific condition 
on the representational accuracy. 

We employ the same basic restructuring algorithms as developed by 
Fritts^ for Cartesian geometry. By conserving the linear momentum, circula¬ 
tion, divergence and the r and z components of the angular momentum on a 
quadrilateral (or triangle) a unique and reversible solution for the new r 
and z components of the triangle velocities is obtained. Conserving the 0 
component of the angular momentum and the generation of circulation gives 
rise to a unique solution for the new angular momentum. These reconnection 
algorithms will be discussed in more detail in a future report. In the next 
section we discuss the equations of motion. 


III. Equations of Motion 

The hydrodynamic equations of motion governing an incompressible, 
inviscid fluid in a cylindrical coordinate system (r, 9, z with corres¬ 
ponding velocity components u, v, w) are 

3u , 3u , 3u v 2 1 3p 

3t 3r 3z r p 3r 


(III-l) 


3v , 3v . 3v . vu . 
—- + U — + w — + — = 0 
3t 3r 3z r 


(III-2) 


3w , 3w 3w 13p 

— + u — + w — =■-, 

3t 3r 3z p 3z 


along with the incompressibility condition 


(II1—3) 
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(III-4) 


i a , . . a A 

7 7F (ru) + a7 w 3 0 


Equation III-2 is just the conservation of angular momentum per unit 
mass. Defining the angular momentum per unit mass L = vr, Equation III-2 is 


3L , m dL 
—+u- VL=—=0 


(III-5) 


We define a pseudo-Cartesian vector space with a gradient operator 

, (III-6) 


rr * - * 3 , « 3 

V' = e t— + e t 
r 3r z 3z 


a vector velocity 


u' = e u + e w 
— r z 


and a position vector 


(II1—7) 


r' = e r + e z 


(III-8) 


Equation III-l and III-3 can then be combined to give 


, d V ' , 

( dc> H. 


1 ~ 
— V p + — e 
p r r 


where 


(III-9) 


<£>’ + V 


Rewriting the r-component of Eq. III-9 we have 

L 2 3£ 
or ~ p 7T = ~ ^ 

This is precisely the equation of motion that would be obtained from a one 
dimensional Lagrangian 

1 1 L" 


V 1. 2 1 L 

=■ 2 pr “ 2 *7* ~ 

with a potential 


(III-10) 


7 


K fTP 










(III-ll) 


V 


p 



L 2 

p r2 


The equations of motion in our pseudo-Cartesian coordinate system become 

(£)’ l -° 

and the angular momentum of a fluid element, per unit mass, remains constant 
as we follow it with its motion; and 


< -{c y - = ■ 0 J ' p " I l2v ' r* • (HI-12) 

The radial and axial motion takes place as if v were absent and, 
instead, a centrifugal force, = - py L 2 V' , were acting in the radial 

direction. Note that the angular momentum is constant in the region over 
which the gradient is applied. It is therefore a triangle quantity — the 
basic fluid element in this code — and each triangle resides in a i/r 2 
potential. The incompressibility condition takes the form 

V' • (ru’) = 0 . (III-13) 

Equations III-ll, 12, 13 now have the same formulation as the 
Cartesian version of the code (SPLISH). The finite-difference algorithms 
and the restructuring algorithms developed for SPLISH can then be applied 
to cylindrical geometry with modifications to include conservation of 
vorticity generation and conservation of angular momentum. The pseudo- 
Cartesian coordinate system will be utilized throughout the remainder of 
this paper. For that reason the primed notation will be dropped in what 
follows. In the next section we discuss the finite-difference algorithms. 
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IV. Finite-Difference Algorithms 


The computer code SPLASH is a 2 1/2-dimensional Lagrangian fluid 
dynamics code for incompressible fluids in cylindrical geometry. This code 
is a direct extension of the philosophy and numerical techniques developed for 
SPLISH, the Cartesian version of the code. As such, most remarks in this 
section apply equally as well to SPLISH. 

The basic equations are. 



(IV-1) 



or explicitly if velocities are considered to be triangle-centered. This 
placement of velocities as cell quantities and pressures at vertices is 
apparently unique to SPLASH and SPLISH and is the direct opposite of the 
usual placement. In what follows the subscript i will denote a vertex- 
centered quantity and j a triangle-centered quantity. In both codes the 
integration of velocities uses a split step algorithm whereby the velocities 
are advanced one half timestep (Eq. IV-4) , the grid is advanced a full 













(IV-6) 


R^ * + 5tu * 


0 j = R[(R°) ,(rJ)] • 0 l 1 


(IV-7) 



^ 2 -^ (7p) rH (7 K 


(IV-8) 


The vertex velocity U_^ appearing in Eq. IV-5 is obtained from the area- 


weighted . from the previous iteration, 

Iu n A. 

■ ~J J 
U n = 1- — - -- 
-i EA. 

J 


(IV-9) 


The advantage of using triangle centered velocities is the ease in concep¬ 
tualizing and expressing conservation laws. Because of the paucity of 
experience in formulating algorithms over a general triangular grid, we 
employed a control volume approach, which uses an integral formulation to 
derive the difference algorithms. Equation IV-7 is the first manifestation 
of this approach. It reflects numerically the fact that the triangle 
velocities must rotate and stretch as the grid rotates and stretches. 

The transformation R is derived by considering the circulation about each 
vertex. The boundaries of a vertex cell are defined by the triangle side 
bisectors as noted in Section II. 

The vertex of Fig. 1 is constructed by summing over all the surrounding 
triangles. Therefore the area of a vertex cell may be defined as 


A =It a - , (IV-10) 

c . 3 j 
J 

where the sum extends over all adjacent triangles. With this definition the 
vertex velocity becomes 
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Since the triangle velocities are constant over the triangle, the circula¬ 
tion taken about the boundary of the vertex cell is straightforward. Circu¬ 
lation is conserved about each of the three triangle vertices by the trans¬ 
formation R of Eq. IV-7. This transformation ensures that the vorticity 
integral calculated about any interior vertex is invariant under the 
advancement of the grid. It is easy to show that the Vp term cannot alter 


the vorticity either since numerically VxVp = 0. Only the (Vp /p.) and 

j J 


(l 2 /2 terms can change vorticity, exactly as dictated by the physics. 

Since the transformation R is time reversible, so are Eqs. IV-4 - IV-3, so 
that the entire algorithm advances vertex positions and velocities reversibly 
while evolving the correct vorticity about every interior vertex. No 
numerical generation of nonphysical vorticity can occur, a rather unique 
feature of both SPLASH and SPLISH among Lagrangian codes. 


The pressure p^ in Eq. IV-8 is derived from the condition that the 


new velocities U . should be divergence free at the new timestep, satisfying 


Eq. IV-2. The pressure Poisson equation is derived from Eq. IV-8 by setting 


(7 - r IJ?)^ = 0, to obtain pressure p^ such that 


Jjt 

2 


r; (Vp, ° 


m (v< r U S - — v 

i 11 j-j'i 2 


r.L? . 


(IV-12) 


where 









Both terms in Eq. IV-12 are simple to evaluate since the divergence is taken 
over triangle centered quantities. The paths are the "surfaces" bounding 
the vertex volume of Fig. 1, where the normal is directed outward from the 
vertex. The Poisson equation (Eq. IV-12) that results from this integra¬ 
tion has two advantages. First it is derived from V 2 <+> = 7*V<(> as in the 
continuum case. Secondly the left-hand side results in the more familiar 
second order accurate templates (such as the five-point formula) for the 
Laplacians for homogeneous fluids and regular mesh geometries. 

In summary the finite difference formulas for SPLASH are derived using 
a control volume approach. Specifications of pressure at vertices leads 
naturally to the choice of positioning velocities at triangles. Although 
pressure gradients are constant over triangles, the resultant algorithms are 
expected to be nearly second order accurate since vertex velocities are 
derived from pressure gradient forces through sums about vertices, which in 
effect centralizes the differences. SPLISH has been tested extensively on 

finite amplitude standing waves and has been shown to be basically second- 

1 12 

order accurate by studying the variation in period with mesh size. * 

V. Couette Flow 

Couette flow refers to the circular flow of a fluid between two rota¬ 
ting coaxial cylinders. This flow is potentially unstable; the instability 
results from a prevailing adverse gradient of angular momentum. 

The stability of an ideal fluid in circulatory motion was first 

13 

investigated by Rayleigh. Simply stated, Rayleigh's criterion says that 
in the absence of viscosity, the necessary and sufficient condition for a 
distribution of angular velocity .'Ur) to be stable is 











~~ (r 2 ft) 2 > 0 

dr 

everywhere in the interval and that the distribution is unstable if (r 2 &) 2 
should decrease anywhere in the interval. 

Note that r 2 ft is the angular momentum, per unit mass, of a fluid element 
about the axis of rotation. We have shown in Section III that the angular 
momentum of an ideal fluid element is a constant of the motion and that the 
motions along the radial and axial direction may be treated as if the 
circulatory motion were absent and instead a centrifugal force [-pL 2 /2V(1/r 2 )] 
were acting in the radial direction. Thus we may associate with each fluid 
element a "potential energy" pL z /2r 2 . This is analogous to the problem of a 
heterogeneous fluid in a field with a potential energy proportional to r“ 2 . 

The equilibrium is stable only if the potential energy is a minimum; i.e., 
the "heavier" fluids are in regions of lower potential energy. This means 

that L 2 must be monotonically increasing outwards. 

3 

Taylor extended this criterion for stability to account for viscosity, 
verified his calculations experimentally and described the secondary flow 
which appears after the onset of the instability. Viscosity tends to 
produce an angular momentum distribution proportional to Ar 2 +B for laminar 
Couette flow, where A and B are two constants related to the angular 
velocities of the inner and outer cylinders. If this distribution is 
unstable, fluid elements with larger angular momenta will move outwards 
inducing a secondary flow. Viscous forces will tend to retard this motion 
but if a viscous drag is not strong enough a redistribution of angular 
momentum will occur. At the same time, the moving solid surfaces will tend 
to re-establish the original distribution of vorticity and a steady 
secondary flow is established. This secondary flow consists of a regular 
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cellular vortex structure in which closed ring vortices alternating in sign 

are wrapped around the axis of symmetry. 

The transiton from steady, laminar Couette flow into fully developed 

Taylor vortex flow is the phenomenon we are simulating here. To obtain a 

quantitative comparison with theory, we compare the growth rates of the 

instability with those obtained from the linear theory developed by 
4 

Chandrasekhar. 

Following Chandrasekhar, we linearize the equations of inviscid motion 
(Eqs. III-l - III-4) by doing a perturbation expansion about the stationary 
flow solution 

u = eu^ 1 ) + e 2 u( 2 ) + ... 

w = ew^ 1 ) + e 2 w( 2 ) + . . . 

v = V°(r) + ev^ 1 ) + e 2 v( 2 ) + . . . 

and p * p°(r) + ep^ 1 ) + e 2 p( 2 ) + . . . 

where V°(r) = rft(r). Assuming all perturbed quantities vary as expi.i(ot-l-kz) ] 

where a is a constant (which can be complex) and k is the wave number of 
the disturbance in the z-direction, we have to first order 


iau^ 1 ) = 2ftv( 1 ) - ~ x~ p^ 1 ^ , (V-l) 

p dr 

iav^ 1 ) + ^ + u^ 1 - 1 = 0 (V-2) 

iavCl) » - — pC1) (V-3) 

P 

+ ikw^ 1 ^ * 0 . (V-4) 

3r r 


combining Eqs. V-l and V-2 and Eqs. V-3 and V-4 and eliminating the 
pressure between the resulting two equations we have (dropping superscripts) 










(V—5) 


where 


[~ 7 ~(ru)] * (k 2 /a 2 )[a 2 - <t(r)]u 
dr r dr 


Mr) - f £ (r 2 G) , 


along with the boundary conditions on the inner (R ) and outer (R 2 ) cylinders 


u(R r ) = u(R 2 ) = 0 


Since the numerical code is Lagrangian, we rewrite Eq. V-5 in terms of 

4 

Lagrangian displacement variable C. 
u = io£ 


v = ia? e - r d? 5 r 


w = iaC 


We also express the angular velocity in terms of its viscid distribution 
il(r) = A + B/r 2 . Equation V-5 then takes the form 


(DD -k 2 )£ * - —* 7 - A(A+B/r 2 )£ 


r a 


(V- 6 ) 


where 


D = j- and D. + - 

dr * dr r 


and the boundary conditions are 


= 0 at r = Rj, R 2 . 


Equation V -6 can be solved explicitly in terms of Airy functions for 
14 

the case of a small gap 

d =* (R t + R 2 ) << 1/2(R x + R 2 ). 


i 









T 


The solutions are complicated in that the wavenumber and growth rate are 

linked parametrically. The coupling terms in the pair of resulting equations 

are determined from a characteristic equation expressed as a ratio of Airy 

functions. We use the growth rates of the most unstable mode as determined 
14 4 

by Reid and Chandrasekhar to compare with the numerical simulation. 


VI. Numerical Results 

The initial grid for the Couette flow simulation is shown in Fig. 3. 
Ri(R 2 ) is the radius of the inner (outer) cylinder. z q (z n ) is the left 
(right) boundary of the computational region. The 8 coordinate is into the 
page. The boundary conditions are 


u(R x ) = u(R 2 ) = 0 

and the system is periodic in the z-direction 


and 


Z o Z N+1 


w(z o ) * W( W 


U(Z 0 } * U( W 


P<z 0 ) 


P<? 

p(z 


N+l ) ’ 
N+l) 


The initial angular velocity is assumed to have the viscid distri¬ 
bution 

»l(r) - A + B/r 2 

where A and B are constants which depend on the angular frequency of the 
inner and outer cylinders. The system is initially in equilibrium with a 
pressure distribution given by 


p°(r,z,t=0) 



Ri 


v 2 (r,z) 

r 


dr + C 
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SPLASH 



Fig. 3 — The initial computational grid for the Taylor-Couette problem. The 
^-direction is into the paper. R 1 (R 2 ) is the radius of the inner (outer) cylin¬ 
der. The r-component of velocity for vertices 1, 2 and 3 are followed in time. 


17 











where C is an arbitrary constant except in the case of a free surface at 
the inner cylinder when it is zero. 

This initial pressure distribution is then perturbed with a 1 % sinu¬ 
soidal perturbation, 

p'(r,z,t=Q) = p°(r,z,t=0) + 0.01p°(r,z,t=Q) sin kz 

where k is the wavenumber of the most unsable mode. The most unstable mode 
is the mode for which the Taylor vortices have a wavelength equal to twice 
the gap width, i.e., the vortices are approximately square in cross-section. 

Two important parameters governing the stability or instability of 
Couette flow are the ratio of radii of the inner and outer cylinders 

n = R^/k^ * 

and the ratio of the angular frequencies of the outer and inner cylinders 

y * 

The Rayleigh criterion for stability can then be written as 

u > n 2 

For all the results presented here, the small gap approximation is valid. 

For Case I we have 

R - 21 cm, R =23 cm, 

1 2 

z„ * 4 cm, k * 2 it/ A ■ ir/2 
N 

and 

* 40tt sec -1 , u ■ 1/2 , 

and this system is unstable since y < n 2 . The initial grid is shown in 
Fig. 3 where we have denoted three vertices which we will follow in time. 








Figure 4 shows the system shortly after one full revolution. As one 
would expect the fluid near the inner cylinder with its larger angular 
momentum is being pushed to a larger radius and the fluid near the outer 
cylinder with its smaller angular momentum is being pushed to a smaller 
radius, i.e., the "heavy" fluid is falling and the "light" fluid is rising. 
Note the large number of vertices that have been added near the boundaries 
to preserve the resolution there. 

The time evolution of the r-component of velocity for the three 
aforementioned vertices is shown in Fig. 5. The growth rates for the three 
vertices are equal for nearly a full revolution at which point vertex 1 
slows down and vertices 2 and 3 speed up. The growth rate in the linear 
regime y c = 213.09 s -i is in very good agreement with the predicted growth 
rate y^ = 215.26 s -1 obtained from linear theory. 

Vertices 2 and 3 are speeding up because they are becoming entrained 
between two very large counter-rotating vortices. Vertex 1 is nearing the 
inner cylinder and being deflected in the z-direction. This is shown 
clearly in Fig. 6 where we have plotted the pathlines of the flow. The 
plus signs are the most recent positions of the vertices and the dots are 
the positions at three previous equal periods of time. The plus signs can 
be regarded as arrow heads. We see the development of two large counter¬ 
rotating Taylor vortex cells. The wavelengths of the vortex cells is 4 cm 
as predicted by the theory. 

Figures 3, 4 and 6 can be directly compared to the experimental work 
of Donnelly and Fultz^ (see Fig. 7). They ejected dye from the inner 
cylinder and took a remarkable series of photographs showing the transition 

from laminar flow to fully developed Taylor vortex flow. Their work is also 

4 

illustrated in Chandrasekhar. 
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PATHL I NES 


rf.lJMlT 



Fig. 6 — Pathlines for Case I. The + signs are the most recent positions of the 
vertices and the dots are the positions at three previous equal periods of time. 
The center vortex has a clockwise rotation. Vertices 1, 2 and 3 are noted. 










The onset of instability with the outer cylinder at rest, n = 0, P c = 4-4'Jl soo: 

(a) laminar flow, P == 4-500 soc ; (b) beginning of radial motion at P — 4-483 sec; (c) Appear¬ 
ance of cells with the outer cylinder at rest; cells at marginal stability, P — P e = 4-466 soc; 

(d) Appearance of cells with cylinders rotating in the samo direction: P — P e = 3-844 sec, 

p = 0-1164. 

Fig. 7 — Photographs showing the transition from laminar Couette flow to fully developed 
Taylor-vortex flow with the outer cylinder at rest. Dye is ejected from the inner cylinder. 
It.J. Donnelly and D. Fultz, Proc. Roy Soc. A 258 , 101 (1960). Used by permission. 


23 




























For Case II, the outer cylinder is at rest, y = 0, and the system 
length has been increased to 12 cm. R^, R^, and k are the same as for 
Case I. Since the wavelength of the perturbation is 4 cm, six Taylor cells 
should develop. This is illustrated in Fig. 8 where we have plotted contours 
of constant circulation. The plus sign indicates flow in the clockwise 
direction. 

The case for which the cylinders are rotating in the opposite direction 
is particularly interesting. Theory predicts that only the inner region of 
fluid should be unstable as it is only in this region that the angular 
momentum is decreasing. For Case III we have 
R 1 = 41 cm, R 2 = 45 cm, 
z^ = 8 cm, k = it 

= 40ns -1 , y = -0.87. 

This choice of y gives zero angular momentum at the center of the gap. The 
constant circulation contours are shown in Fig. 9. As predicted, four 
counter-rotating vortex cells appear in the inner half of the fluid while 
the outer half of the fluid remains stable. This figure can also be 
directly compared with the photographo of the experimental work of Donnelly 
and Fultz 1 '* (Fig. 10). 

The numerical results for all the cases are summarized in table form 
in Fig. 11. The computational growth rates are in very good agreement with 
the theoretical growth rates with errors on the order of 1%. Note that for 
u = 1 the fluid is stable and y is the oscillation frequency of the 
perturbed velocity. 
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NUMERICAL RESULTS 
COUETTE FLOW 

V = 0.91 


H- 

£l t (rad/i) 

k 

--- 
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2 0 ir 

7T / 3 
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Fig. 11 — Comparison of computational growth rates with growth rates obtained 
from linear theory for various initial conditions 
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VII. Summary and Conclusions 

We have developed a 2 1/2 dimensional, Lagrangian, hydrodynamic model 
designed for Che simulation of transient rotational flow phenomena. The 
model uses as a finite difference mesh a general connectivity triangular 
grid. The advantages of this model are numerous: 

1) Complicated geometries, interfaces and free surfaces can be 
treated with a minimum of difficulty. 

2) The resolution across the mesh is highly variable. 

3) The mesh can be restructured to preserve the numerical accuracy 
of the simulation. 

4) No numerical vorticity is generated. 

5) The gradient operator, as a triangle function, is exact in the 
linear approximation. 

6) The integral operator, as a vertex function, is exact in the linear 
approximation. 

We have applied this model to the study of the transition of laminar 
Couette flow to Taylor vortex flow with a high degree of success. The 
computational growth rates are in excellent agreement with the theory. 

For the results presented here the code has not been run long enough 
to achieve steady-state Taylor-vortex flow. However, with respect to 
the transition to Taylor-vortex flow we find that a vortex signature appears 
very early in the run when the system is still in the linear regime. By 
this we mean that contours of constant circulation show uniformly spaced 
vortices when the perturbed velocities are on the order of 10~ 6 cm/sec. 

[V. 0(10 3 cm/sec)]. These vortices then increase in strength but maintain 


their shape and spacing. 








In these calculations we have perturbed the system at only one wave¬ 
length corresponding to the predicted wavenumber for the steady-state Taylor 
vortices. Although the wavenumber is therefore not expected to change, it 
is surprising to find that the flow throughout the cylinder gap is estab¬ 
lished during the linear regime in exactly the nonlinear flow pattern. As 
evidenced by the case of the counter-rotating cylinders, this flow is not 
evident a priori. The transition to Taylor-vortex flow from laminar 
Couette flow is strickly speaking nonexistent— the linear flow is only 
strengthened. Whether any consequences of this simple transitioning are 
evident in the more complicated cases of viscid flow a d the perturbations of 
many wavenumbers will be investigated in future calculations. 

Two code modifications are necessary for these calculations to be 
compared with experiments; the addition of viscosity (see Appendix A) and 
the improvement of some grid restructuring algorithms. In order to 
preserve numerical accuracy when the grid becomes highly distorted vertex 
additions and deletions must be made. We have now developed vertex 
addition/deletion algorithms which conserve divergence, curl, linear 
momentum, angular momentum and vorticity generation. These algorithms are 
now being incorporated into the code to allow the calculation to proceed 
further into the fully nonlinear regime while preserving second-order 
accuracy. These vertex addition/delection algorithms do not affect the 
results presented here. 

The inclusion of viscosity and the new restructuring algorithms will 
enable a direct comparison between the computational results and the 
experimental results to be made. This would entail calculating the torques 
required to maintain the cylinders in motion and the critical Taylor number, 
the ratio of the destabilizing centrifugal force and the opposed radial 
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Appendix A 


Consider the viscid equations of motion written in our pseudo- 
Cartesian format 

4“ . - I Xl u f 2 J./3u) + J_ /3u + Sw) + 2/3u _ u\] ( 

dt p ar + r + o[ 2 3rl3r/ 3z Uz + 3r/ + r(3r r/J* (A 1} 


and 


dw 

dt 


dL 

dt 


= _Ii2 + M./ 2 — [—1 + - — Fr/— + —\11 (A-2) 

p 32 p 1 32 1.3zj r 3r L r \3z 3r/Jj ' ' 


f 3 2 v , 3 / 3v v\ 2 / 3v v\1 

rp L3z 2 3r \ 3r - r) r \ 3r r/J 


The coefficient of viscosity p is considered a constant. 

A-2 can be combined to give 

du - 

5J---7P- 1/21^-£±7x0*) , 

where uj_ is the theta component of vorticity 


(A-3) 

Equations A-l and 

(A-4) 


uj = u) a e = Vxu 
o y 

u) is a vertex function which is easily determined by calculating the circula- 

0 

tion about a vertex 


u», = ~ § u’dl 

8 i A i “ ~ 


Equation A-3 can be cast into the form of a diffusion equation for the 
angular momentum. 

^-ryV-^VL) . (A-5) 

The finite differencing of this equation is discussed in detail in Ref. 10. 
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